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First applications of the new algorithm simulating dynamical fermions are reported. The method reproduces 
previous results obtained with different techniques. 



1. DYNAMICAL FERMIONS AND 
LUSCHER REPRESENTATION 

The problem of dynamical fermions, even 
though solved in principle, still provides the con- 
siderable challenge in practice. Existing exact al- 
gorithms require huge computing resources due 
to the strong correlations between generated con- 
figurations Therefore the recent proposal of 
Liischer has attracted a lot of interest J2j , (D . To 
remind, the main difficulty with simulating the 
theory with dynamical fermions consists of the 
nonlocality of the fermionic determinant 



[dH>dH>] exp - det(M), 



(1) 



which depends functionally on the gauge field in 
the case of QCD for example. The local and 
positive representation for the positive powers 
of det(A4) was not found up to date. Liischer 
method exploits the well known fact that the in- 
verse of the determinant has such a represente- 
tions in terms of the standard bosonic fields. He 
therefore proposes to use the polynomial approx- 
imation of the reciprocal function 



- = lim Pn(s). 
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(2) 



In the simplest case P/v(s) can be taken as the 
geometrical series (in 1 — s), in practice better 
choices are known. Decomposing Pn into the 
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product form and using the symmetries of roots 
of real polynomials Liischer proves that the par- 
tition function of the full QCD with two flavours 
of dynamical quarks has positive and local repre- 
sentation in terms of the original gauge fields and 
N bosonic fields associated with zeros of P/v(s). 

It remains now to be tested if this elegant trick 
is also advantageous in practice. In particular 
the crucial question is: how do the parameters 
of this representation scale with the volume and 
with the bare parameters of the original theory in 
its critical region? All these important and crit- 
ical problems are being now thouroghly checked 
1- 

In this contribution I report on the first re- 
sults of applying the Liischer 's method to the 
(Euclidean) three-dimensional Hubbard model, 
which is the convenient and nontrivial testing 
ground of the new techniques dealing with dy- 
namical fermions. 



2. THE HUBBARD MODEL AND ITS 
BOSONIC REPRESENTATION 

The Hubbard model is a many-body theory of 
nonrelativistic electrons with the nonlinear point- 
like coupling representing the effective Coulomb 
interactions §,(6). Its particular Euclidean for- 
mulation contains also the continuous Hubbard- 
Stratonovich field which couples to the electrons, 
hence may be regarded as analogous to the Yang- 
Mills fields of QCD. Even though the nonrela- 
tivistic, the model contains all the numerical dif- 
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Acuities caused by the fermionic statistcs, i.e. by 
the Pauli exclusion principle. Being simpler than 
QCD it still describes interesting physics of a non- 
trivial many-body system. 

Hubbard Hamiltonian may be written as 

U 



h = -k a L a i* - j ~ n ^"> 2 ( 3 ) 

<ij>a i 

(4) 



where 



- „t 



ailn-Qio-i an d a-icr denotes the creation 



operator of an electron at the lattice site i and 
with the spin a =|, !• K, U and \i are the hopping 
parameter, coupling constant and the chemical 
potential respectively. Standard transfer matrix 
formalism gives the following Euclidean represen- 
tation for the partition function |7j] 

Z= f [dA]exp (-^2A%/2) detM+detM-,(5) 

where the continuous Hubbard-Stratonovich field 
An was introduced to decouple the quartic in- 
teraction term. The discrete index t — l,...,Nt 
labels the time slices. Fermionic matrices M± 
define the bilinear fermionic actions 



<ij>t 



+ J2^u (exp 

it \ 



I^A it - { U±,)L 
N t V w JV t 
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At fi = 0, which corresponds to the half-filling, 



detM+dctM- = dctM = det MM 



(6) 



because the determinants are real. We can there- 
fore apply the Liischer trick and write the poly- 
nomial approximation for the inverse 



1 



2/V 



oto - P2n{Q ] Q) = JJ(Q f Q - z k ) 

N 



(7) 
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since the roots z k — a k + iftk of the polyno- 
mial P2n{z) exist in the complex conjugate pairs. 



Q f Q = M^M/Xmax with \ max being the largest 
eigenvalue of M'M. 

We therefore have the local and positive repre- 
sentation for the partition function of the Hub- 
bard model. 



Z ~ J [dAd<j>] exp (- J A 2 (x)/2 d 3 a 

ex p(" [{Q^Q- ot-kf + Pi) 

\ J fc=i 



l ) k d 3 x 



(8) 



(9) 



where we have used continuous notation for sim- 
plicity. Hence the original system of nonlinearly 
interacting fermions was replaced by TV complex 
bosonic fields coupled to the single real scalar 
field A(x). This can be simulated with the stan- 
dard local Monte Carlo techniques. Compared 
to the original Liischer mapping we have half as 
many fields 4> k (x), however our bosonic action is 
more complicated since Q} Q contains also next- 
to-nearest- neighbour interactions. 

2.1. IMPLEMENTATION AND 
RESULTS 

Convergence of the polynomial approximation, 
Eq.(^), cannot be uniform in s in the whole in- 
terval (0,1). Hence Liischer introduces the cut-off 
e > such that the series (||) is uniformly conver- 
gent for e < s < 0. This parameter enters directly 
into the construction of the optimal polynomials 
Pn ■ The value of e is crucial for the practical ap- 
plicability of the algorithm since it controls the 
number of bosonic fields required to reach the pre- 
scribed accuracy. Ideally e should be smaller that 
the smallest eigenvalue of the normalized matrix 
Q^Q- Otherwise corrections for a few eigenvalues 
lower than e could be implemented || . 

Our preliminary simulations were done on the 
5 3 lattice at K = 1 and U = 1. These param- 
eters were chosen to allow comparison with the 
earlier Creutz results obtained with different al- 
gorithm . One Monte Carlo sweep consisted of 
the heatbath generation of all fields <j) k and the 
Metropolis update of the Hubbard-Stratonovich 
field A(x). Fig. 1 shows the distribution of the 
minimal eigenvalue of Q) Q. On this basis we have 
chosen e = 0.001. The number of auxiliary fields 
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Figure 1. Distribution of the smallest eigenvalue 
of the matrix Q) Q found in the simulation de- 
scribed in text. 



4>k was fixed to N = 100. This guaranteed that 
the relative error introduced by finite N in Eq. ^ 
was smaller than 10 ~ 4 for e < s. Instead of X max 
needed to define the normalized matrix Q we have 
used the upper bound A which can be readily de- 
rived 

X max < 2[X max (A^A) + max(D) 2 } = A, (10) 

where D is the diagonal part of M. — A + D. 
Hubbard- Stratonovich field was cut-off by A(x) < 
4 which would correspond to the four standard 
deviations if the determinants were neglected. 
This bound was never violated in practice. Fig. 2 
shows the history of the density of the electrons 
n f and of the local number of pairs ri](x)ni{x). 
Both observables stabilize relatively quickly and 
the equilibrium averages agree with the results 
quoted by Creutz |?J for slightly bigger lattice. 
For the 6 2 x 8 lattice we obtain < n T >= 0.47(3) 
and < ntni >= 0.220(3) while the correspond- 
ing numbers read from the Creutz plots are 0.48 
and 0.22 with the read-out error ~ 0.01. The 
run summarized in Fig. 2 took about 10 hrs of the 
HP-720. 

We have found however that the configurations 
generated by this algorithm are strongly corre- 
lated. The autocorrelation time measured on 
both densities is approximately 100 sweeps. Sim- 
ilar phenomenon was also reported at this Con- 
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Figure 2. The history of run of total length 3000 
events. The larger values are for density of elec- 
trons, the smaller for density of pairs. 



ference by Jegerlehner while testing QCD imple- 
mentation M. 

In conclusion, the Hubbard Model is in the 
Luscher class, i.e. it can be mapped onto a system 
of bosons with local interactions. The half- filed 
case has a positive Boltzman factor which admits 
standard Monte Carlo approach. Further studies 
are needed to investigate whether the large auto- 
correlation times can be eliminated. 
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